#include "BicubicUniformPeriodicBSplinePatch.h"

using namespace cagd;

BicubicUniformPeriodicBSplinePatch::BicubicUniformPeriodicBSplinePatch(
        GLdouble u_min, GLdouble u_max, GLdouble v_min, GLdouble v_max,
        GLuint row_count, GLuint column_count, GLboolean is_closed):
        AbstractTensorProductSurface(u_min, u_max, v_min, v_max,
                                     row_count, column_count, is_closed)
{}

BicubicUniformPeriodicBSplinePatch::~BicubicUniformPeriodicBSplinePatch()
{}

GLboolean BicubicUniformPeriodicBSplinePatch::UBlendingFunctionValues(
        GLdouble u_knot,
        RowMatrix<GLdouble> &blending_values) const
{
    if (u_knot < 0.0 || u_knot > 1.0)
        return false;

    blending_values.ResizeColumns(4);

    GLdouble u = u_knot, u2 = u * u, u3 = u2 * u, w = 1.0 - u, w2 = w * w, w3 = w2 * w;

    blending_values(0) = w3 / 6.0;
    blending_values(1) = 0.5 * (u * w2 + w) + 0.16666666666666666666666666666667;
    blending_values(2) = 0.5 * (u2 * w + u) + 0.16666666666666666666666666666667;
    blending_values(3) = u3 / 6.0;

    return true;
}

GLboolean BicubicUniformPeriodicBSplinePatch::VBlendingFunctionValues(
        GLdouble v_knot,
        RowMatrix<GLdouble> &blending_values) const
{
    if (v_knot < 0.0 || v_knot > 1.0)
        return false;

    blending_values.ResizeColumns(4);

    GLdouble v = v_knot, v2 = v * v, v3 = v2 * v, w = 1.0 - v, w2 = w * w, w3 = w2 * w;

    blending_values(0) = w3 / 6.0;
    blending_values(1) = 0.5 * (v * w2 + w) + 0.16666666666666666666666666666667;
    blending_values(2) = 0.5 * (v2 * w + v) + 0.16666666666666666666666666666667;
    blending_values(3) = v3 / 6.0;

    return true;
}

GLboolean BicubicUniformPeriodicBSplinePatch::CalculateAllInformation(
        GLdouble u,
        GLdouble v,
        SurfaceInformation &pd) const
{
    if (u < 0.0 || u > 1.0 || v < 0.0 || v > 1.0)
        return false;

    RowMatrix<GLdouble> u_blending_values(4), d1_u_blending_values(4);

    GLdouble u2 = u * u, u3 = u2 * u, wu = 1.0 - u, wu2 = wu * wu, wu3 = wu2 * wu;

    u_blending_values(0) = wu3 / 6.0;
    u_blending_values(1) = 0.5 * (u * wu2 + wu) + 0.16666666666666666666666666666667;
    u_blending_values(2) = 0.5 * (u2 * wu + u) + 0.16666666666666666666666666666667;
    u_blending_values(3) = u3 / 6.0;

    d1_u_blending_values(0) = -0.5 * wu2;
    d1_u_blending_values(1) = 0.5 * (wu2 - 1) - u * wu;
    d1_u_blending_values(2) = wu - 0.5 * (u - 1.0);
    d1_u_blending_values(3) = 0.5 * u2;

    RowMatrix<GLdouble> v_blending_values(4), d1_v_blending_values(4);

    GLdouble v2 = v * v, v3 = v2 * v, wv = 1.0 - v, wv2 = wv * wv, wv3 = wv2 * wv;

    v_blending_values(0) = wv3 / 6.0;
    v_blending_values(1) = 0.5 * (v * wv2 + wv) + 0.16666666666666666666666666666667;
    v_blending_values(2) = 0.5 * (v2 * wv + v) + 0.16666666666666666666666666666667;
    v_blending_values(3) = v3 / 6.0;

    d1_v_blending_values(0) = -0.5 * wv2;
    d1_v_blending_values(1) = 0.5 * (wv2 - 1) - v * wv;
    d1_v_blending_values(2) = wv - 0.5 * (v - 1.0);
    d1_v_blending_values(3) = 0.5 * v2;

    pd.LoadNullVectors();
    for (GLuint row = 0; row < 4; row++)
    {
        DCoordinate aux_d0_v, aux_d1_v;
        for (GLuint column = 0; column < 4; column++)
        {
            aux_d0_v += _control_net(row, column) * v_blending_values(column);
            aux_d1_v += _control_net(row, column) * d1_v_blending_values(column);
        }
        pd.point    += aux_d0_v * u_blending_values(row);
        pd.diff1u   += aux_d0_v * d1_u_blending_values(row);
        pd.diff1v   += aux_d1_v * u_blending_values(row);
    }

    return true;
}
